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The basic purpose of the paper is to draw the attention of researchers to new possibilities of 
differentiation of similar signals having different nature. One of examples of such kind of signals is 
presented by seismograms containing recordings of earthquakes (EQ's) and technogenic explosions 
(TE's). EQ's are among the most dramatic phenomena in nature. We propose here a discrete 
stochastic model for possible solution of a problem of strong EQ's forecasting and differentiation of 
TE's from the weak EQ's. Theoretical analysis is performed by two independent methods: with the 
use of statistical theory of discrete non-Markov stochastic processes (Phys. Rev. E 62,6178 (2000)) 
and the local Hurst exponent. The following earth states have been considered among them: before 
(lb) and during (I) strong EQ, during weak EQ (II) and during TE (III), and in a calm state of Earth 
core (IV). The estimation of states I, II and III has been made on the particular examples of Turkey 
(1999) EQ 's, state IV has been taken as an example of underground TE. Time recordings of seismic 
signals of the first four dynamic orthogonal collective variables, six various plane of phase portrait 
of four dimensional phase space of orthogonal variables and the local Hurst exponent have been 
calculated for the dynamic analysis of states of systems I-IV. The analysis of statistical properties 
of seismic time series I -IV has been realized with the help of a set of discrete time-dependent 
functions (time correlation function and first three memory functions), their power spectra and 
first three points in statistical spectrum of non-Markovity. In all systems studied we have found out 
bizarre combination of the following spectral characteristics: the fractal frequency spectra adjustable 
by phenomena of usual and restricted self-organized criticality, spectra of a white and color noises 
and unusual alternation of Markov and non-Markov effects of long range memory, detected earlier 
(J. Phys. A 27, 5363(1994)) only for hydrodynamic systems. Our research demonstrates that discrete 
non-Markov stochastic processes and long- range memory effects play a crucial role in the behavior of 
seismic systems I -IV. The approaches, permitting to obtain an algorithm of strong EQ's forecasting 
and to differentiate TE's from weak EQ's, have been developed. 

PACS numbers:0.5.65.+b, 02.50.-y, 05.40.-a, 05.65.+b, 45.65.+k 



I. INTRODUCTION 

The earthquakes are among the most mysterious and dramatic phenomena, occurring in nature. As a result of 
sets and breakups of terrestrial cortex or higher part of mantle over hundreds of thousands underground pushes and 
fluctuations of the Earth surface occur annually. They propagate over long distances in the form of elastic seismic 
waves. Nearly thousands of them are registered by people. Annually nearly hundred earthquakes cause catastrophic 
consequences: they affect big communities of people and lead to great economical losses. 

For the study of basic mechanisms underlying its nature modern numerical and statistical methods are used now in 
modeling and understanding of the EQ phenomenon. In papers |l[] , Q the modified renormalization group theory with 
complex critical exponents has been studied for implications of EQ's predictions. Long-periodic corrections found fit 
well the experimental data. Then universal long-periodic corrections based on the modified renormalization group 
theory have been used successfully Q for possible predictions of the failure stress phenomenon foregoing an EQ. The 
failure stress data are in a good reliability with acoustic emission measurements. In paper (Q it has been shown that 
the log-periodic corrections are of general nature, they are related to the discrete scale invariance and complex fractal 
dimension. This idea has been checked in ||, || for diffusion- limited-aggregate clusters. The paradox of the expected 
time until the next EQ with an attempt of finding of acceptable distribution is discussed in fin. New explanation 
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of Guttenberg-Richter power law related to the roughness of the fractured solid surfaces has been outlined in [||. 
Recent achievements and progress in understanding of the complex EQ phenomena from different points of view are 
discussed in the recent review New numerical methods like wavelets and multi-scale singular-spectrum analysis 
in treatment of seismic data are considered in jio) . 

All these previous methods are developed for understanding the statistical and non- stationary properties of EQ's 
and TE's. But in this paper we would like to demonstrate some possibilities related initially to differentiation of EQ's 
from TE's. This problem has not only scientific significance related with recognition of similar signals having physical 
origin. In recent time it was related also with some political problems associated with testing of nuclear explosions 
also. 

Seismic data are an object of careful analysis and numerous methods of their treatment are used especially for 
forecasting of EQ's with strong magnitudes. In spite of wide application of approaches based on nonlinear dynamics 
methods, the Fourier and wavelet transformations etc., we have essential limitations, which narrow down the range 
of applicability of the results obtained. One of the main difficulty is that the discrete character of the seismic signals 
registration is not taken into account. Another factor, which should be taken into account is related to the influence 
of local time effects. Alongside with of the discreteness and the local behavior of the seismic signals considered here 
exist the third peculiarity, viz., the influence of long-range memory effects. 

In this paper, we present one of the possible solutions to forecasting strong EQ's and differentiating TE's from 
weak EQ's. In this presentation we consider of three important factors for seismic signals registered in the form of 
seismograms : discreteness, long- range memory and local time behavior. Two new methods are used to analyse 
these three factors. The first one is based on seismograms considered in the form of a discrete non-Markov statistical 
process along with analysis of corresponding phase portraits, memory functions and the non-Markovity parameters. 
The second method is based on the generalized conception of the Hurst exponent. These methods have been used 
for careful analysis of seismic data and to differentiate EQ's from TE's . The results obtained with the use of these 
methods are useful in recognition of specific features of EQ's and TE's and can be used for strong EQ's forecasting. 

The paper is organized as follows. In Section 2 we describe in brief the stochastic dynamic of time correlation in 
complex systems containing seismic signals by the discrete non-Markov kinetic equations. The local fractal dimension 
and the corresponding Hurst exponent are defined in the Section 3. The real data treatment with the use of non- 
Markov conceptions has been realized in the Section 4. The Section 5 contains some results obtained by the local 
Hurst exponent method. The basic conclusions are discussed in the final Section 6. 



II. THE KINETIC DESCRIPTION OF DISCRETE NON-MARKOV RANDOM PROCESSES 

In recent paper [[fl| the statistical theory of discrete non-Markov random processes has been developed. The basic 
elements, which are necessary for understanding of other sections, are presented, in brief, here. In accordance with 
the Refs. [0- [[l3| the fluctuations of random variable Sxj = Sx(T + jr), j = 0, 1, • • • , N — 1 of a complex system can 
be represented as /c-component state vector 

A£(0) = (5x 0) Sx!,5x 2 , • • • , Sxk-i) = (Sx{T), 6x(T + r), • • • , 6x(T + (k - l)r). (I) 

Here r is a finite discretization time, Sxj and < x > define fluctuations and mean value correspondingly, T is the 
beginning of the time series. They are defined by conventional relationships 

1 

5xj = Xj- <x>, < x >= — x(T + jr). (2) 

The set of state vectors forms a finite-dimensional Euclidean space, where the scalar product of two vectors can be 
defined as 

fe-i 

<A.B>=J2^B j . (3) 

j=o 

The time dependence of the vector A can be defined as result of discrete m-step shift 

A™ +i .(t) = {5x m , 5x m +i,Sx m+ 2, • • ■ , Sx m+ k-i} 
= {Sx(T + mr),Sx(T + (m + l)r), Sx(T + (m + 2)r), • • • , Sx(T + (m + k- l)r}, (4) 



where t — mr and r is a finite time step. Statistical parameters (absolute and relative variances) can be expressed 
by means of the scalar product of two vectors as following 

<ji = ±<A%-A%>=N^{A%y, 

<r A ■ A > 

^2 _ 2: N A N ^ 

N<X> 2 

We define the evolution operator for the description of evolution of the variables dxj as following 

Sxj+i (T + (j + l)r) = U(T + (j + 1)t,T + jr)5 Xj (T + jt) = U(t)5x 3 . (5) 
One can write formally the discrete equation of motion by the use of operator U(t) in the form 

Ax(t) x(t + T)-x(i) 1 



At 



-{U(t + T,t) - l}x(t). (6) 



The normalized time correlation function (TCF) can be represented by Eqns. ([!]) and (|^) (where t = mr is current 
discrete time) as follows 

m _ < A° fc ■ AZ +k > _ <Ag(0)-A^ +fc (t)> 
a[) <Al-A Q k > <A£(0) 2 > ' () 

From the last equation (^) one can see that TCF a(t) is obtained by projection of the final state vector A^ +fc (t) (^) 
on the initial state vector A°(0). Because of this property one can write the projection operator in the linear space 
of state vectors 

< AOfn'lA™ 1 (t) > 

UAZ +k (t) = Ag(0) ^Agfft = A °(°) fl W- ^ 
The projection operator II has the following properties 

n= ^^)f>^ ' n2 = n ' p = 1 " n ' p2 = p ' np = o,pn = o. (9) 

The projection operators II and P are idempotent and mutually complementary. Projector II projects on the direction 
of initial state vector A°(0), while the projector P projects all vectors on the direction which is orthogonal to the 
previous one. Let us apply the projection technique in the state vectors space for deduction of the discrete finite- 
difference Liouville's equation 



A 
At 

L(t,r) = (ir)- 1 {C/(t + r,t)-l}. 



— A£ +k (t)=iL{t,T)A% +k {t), (10) 



The first expression defines the Liouville's quasioperator L and the second expression defines the evolution operator 
U(t). Transferring from vectors A™ +fe to a scalar value of the TCF a(t) by means of suitable projection procedure 
one can obtain the closed finite-difference equation for the initial TCF 

-^l = Aio(t) - rA x J2 Mt{jr)a{t-3T). (11) 

3=0 

Here Ai is the relaxation parameter while the frequency Ai defines the eigen spectrum of the Liouville's quasioperator 
L in the following way 

_ <A°(0)LA°(0)> _ < A°£ 12 L 21 A°(0) > 
Xl - 1 <|Ag(0)P> ' <|A°(0)|*> ' (12j 

The standard Liouville's equation is obtained easily from Eqs. (6), (10) by means of the limit r — » 0. In this case 
the Liouville's quasioperator L is reduced to a classical or quantum Liouvillian and is defined correspondingly by the 



classical or quantum Hamiltonian of the system considered. The given approach is true for non-Hamiltonian systems 
of arbitrary nature when the Hamiltonian cannot be written together with conventional equations of motion. The 
function M\(Jt) in the right hand side of Eqn. (pd|) is the first order memory function 



(13) 



<AUO)L 12{ l^rL 22 yL 21 Al(0)> 
<A°(0)Z 12 Z 21 A°(0)> 



Here we use the following notation for the matrix elements of the splittable Liouvillian quasioperator Li 



ILil!,-, i = 



1,2, ni = II, II 2 = P, L n = nLn, L 12 = ULP, L 21 = PLU, L 22 = PLP. Equation Q can be considered as 
the first equation of the finite-difference kinetic equations chain with memory for the discrete TCF a(t). In paper 
jllj it has been demonstrated that using Gramm-Schmidt orthogonalization procedure one can define the dynamic 
orthogonal variables W„ (t) by means of the following recurrence relationships 



W o = A°(0), Wi ={i£-Ai}W , W n = {ii-A„_ 1 }W n _ 1 + A n _ 1 W n _ 2) n > 1. 
Here we introduce the fundamental eigen values A n and relaxation A n parameters as follows 

_ . < W n XW n > _ _(W„_l(iL-A„+l)W n ) 



< IWJ 2 > 



(IWn-il 2 ) 



(14) 



(15) 



Parameters X n are very similar to the Lyapunov's exponents. Arbitrary orthogonal variables W„ can be expressed 
directly via the initial variable Wp = A"(0) by means of (|l4|) in the generalized form 



(iL-Ar) 



A 



1/2 



k\' 2 (iL - A 2 ) 



A 



1/2 
2 





A l/2 

(iL - A 3 ) 








(iL - A n _i) 



W . 



(16) 



The physical sense of the new variables W„ can be interpreted as follows. For example, the local density of fluctuations 
in the physics of continuous media can be identified with the initial variable Wo. In this case the fluctuations of 
the local current density, local energy density and local energy current density can be associated with the dynamic 
variables W„ with numbers n — 1,2,3, correspondingly. 

One can relate to the set of projection operators n„ to the set of orthogonal variables (|l4]). The last ones project 
an arbitrary dynamic variable (viz., a state vector) Y on the corresponding initial state vector W„ 



n„ 



|Wn ><W.*| 
< IWJ 2 > 



rr 2 — TT P — 1 — TT P 2 — P n P — f) 

HnH m — ^n,mH n , PnPm — ^n.m-^ri; Pn^-n 0. 



(17) 



Acting successively by projection operators n„ and P n on the finite-difference equations ( |To| ) for the normalized 
discrete memory functions 



M n (t) 



< W„[l + iTL 2 2 ) ] m W„ > 



< |W„(0)| 2 > 

one can obtain a chain of the coupled non-Markov finite-difference kinetic equations of the following type 

AM n (t) 



At 



m—1 

X n+1 M n (t) - rA„ +1 V M n+1 



3=0 



(jT)M n (t-jr) , 



(18) 



(19) 



Here X n is the eigen value spectrum of Liouville's operator iL, while A„ is the general relaxation parameters, 

_ . < W ? * t £W» > _ _<W„_i(ii-A„+i)W„) 



< iw n i 2 > 



(|W„_i| 2 ) 



which were defined before by relationships (15). One can consider the set of the functions M n (t) together with the 
initial TCF (n = 0) 



M (t)=a(t) = | A o( Q )|2 > ' = TYIT : 

as functions characterizing the statistical memory of the complex system with discrete time. The initial TCF a(t) and 
the set M n {t) of discrete memory functions appearing from Eqns.(|19|) are playing an important role for the description 
of non-Markov and long-range memory effects. Now it is convenient to rewrite the set of Eqns. (|l^) as the chain of 
the couplednon-Markov discrete equations for initial discrete TCF a(t) (t = wit) and represent them in the form 



Aa(t) 
At 

AMi(t) 
At 

AMj(t) 
At 



Aio(t) - rAi Mi(jT)a(t - jr), 

3=0 

771 — 1 

A 2 Mi(t) - rA 2 M 2 (jT)Mi(t - jr), 
j=o 

771 — 1 

A 3 M 2 (t) - rA 3 J2 M 3 ^t)M 2 (t - jr) . (20) 



3=0 



The kinetic finite-difference Eqns. ( |l9| ) and (^0|) are analogous to the well-known chain of kinetic equations of 
the Zwanzig'-Mori's type. These equations are playing a fundamental role in the modern statistical mechanics of 



nonequilibrium phenomena with continuous time. One can consider the kinetic equations (20) as a discrete-difference 
analogy of hydrodynamic equations for physical phenomena with discrete time. On the basis of the initial set of the 
experimental data one can find the set of orthogonal variables W n in the following way 

W = A°, Wx = ( A - Xi) A°, 



W 2 



At 

A - A 2 ) Wi + A x A° fe = { (£\ - A( Al + A 2 ) + AiA 2 + Ar}A° 



w 3 = (A_ A3 )w 2+ A 2 (A_ Al ) A o. (21) 

It seems to us that one could suggest more physical interpretation of different terms in the right side of the three 
Eqs. (21). For example, term can be associated with dissipation, term is similar to inertia and term AiA(t) 

is related to restoring force. Then the third finite-difference derivative is associated with the finite-difference 
form of the Abraham-Lorenz force corresponding to dissipation feedback due to radiative losses as seen from recent 
experimental evidence in frictional systems |llf| . 

In concrete applications it is necessary to take into account that the dimension of new state vectors W„ is gradually 
decreasing with the increase of the number n. If the initial vector A° has dimension the k then the vectors Wi, W 2 
and W 3 will have the dimensions k — 1, k — 2 and k — 3, correspondingly. 

Solving the chain of Eqns. ( |l9[ ) under the assumption that all A s = one can find recurrence formulae for memory 
functions of arbitrary order in the following form 

Tit — 1 

M s {mr) = - M s (jT)M s ^((m + 1 - j)r) + eJ 1 {M.- 1 {(m + l)r) - M a _ x ((m + 2)r)}, (22) 

3=0 

e s =t 2 A s , s = 1,2,3,- ... 

By analogy with fll|| it is convenient to define the generalized non-Markov parameter for frequency-dependent case 
as follows 

where z = 1,2,,, and (J-i(uj) is the power spectrum of the i-th memory function. It is convenient to use this parameter 
for quantitative description of long-range memory effects in the system considered together with memory functions 
defined above. 

The set of new parameters describes the discrete structure of the system considered and allows to extract additional 
information related to non-Markov properties of the complex (non-Hamilton) systems . 



III. LOCAL HURST DIMENSION ANALYSIS OF SEISMIC DATA 



The Hurst exponent. Typical seismic data are seismic waves registration written in the form of vibrations 
of the Earth surface. Many observations as seismograms lead to random series registrations: technogenic noises, 
gravimetrical, economical, meteorological and other data. Some properties of such random series can be characterized 
by the Hurst exponent H [^9|, (^0|. Let & define the i-th value of the observable variable, < £ T > define its mean value 
on the segment containing r registered points. For the cumulative average value we have X(t, t) = 53i=i(& — < 6 >r)- 
The range R for the given sampling of the random series considered is defined as follows 

R(t) = maxX(t, r) - minX(t, r), (24) 

at 1 < t < r, where t is discrete time accepting integer values and r is a length of the time sampling considered. 
Normalizing the range R on the standard deviation S for the chosen sampling ^ 

S(T)=U y £W)-(t)r}* S J (25) 

and analyzing the variations of the normalized range Hurst p9| , pofl obtained the following empirical relationship: 

S(t) 

where R is the range, S is the normalized variance and H is the so-called Hurst exponent for the sampling of the given 
length £. The value H = 0.5 corresponds to the normal distribution sampling, other values correspond to the various 
degrees of correlations, which can be interpreted in the terms of persistent coefficient. One can use the normalized 
range method for the definition of the Hurst exponent, but it works well for large samplings containing 1000-10000 
registered points. 

The calculation of the Hurst exponent for seismic data. One can obtain easily the Hurst exponent for long 
(1000 — 10000 registered points) samplings |Q by means of the method of the normalized range (R/S analysis). The 
Hurst exponent restoration accuracy calculated on the model data is located in the interval (0.1-1%). For example, 
if the model Hurst exponent was chosen as 0.7 then in the result of the R/S analysis the restored value is equaled to 
0.69 with the changeable third decimal point. The calculated Hurst exponent for the initial seismic noise without an 
"event" (earthquake/explosion) accepted the values 0.96 — 0.98. The obtained values show the high level of persistency 
and correlation. However, these values can be referred to the whole series and cannot reflect the peculiarities of the 
event. In other words, the values of the Hurst exponents calculated for the whole series cannot provide information 
about possible EQ's/TE's, which can be characterized by other values of persistency. In this situation it is necessary 
to generalize this parameter and define the notion of local Hurst exponent. 

The local Hurst exponent. The generalized (local) Hurst exponent can be a sensitive indicator, which gives 
additional information about the regular component in the sampling considered. But the reason for changing of the 
Hurst exponent H is not only the presence of the signal in the sampling considered but slow (for natural processes) 
variations of the correlated noise itself. 

If one considers random scries for a relatively long time it is logically appropriate to cut the series into short segments 
and calculate the Hurst exponent H for each of them. In such a manner, one can detect the variations of H on time or 
in some spatial coordinates. It is better to use the shortest intervals possible for calculating the local exponent H(t). 
Suf ficient number of registered points can serve as a criterion for choosing the minimal interval for such kind of 
statistical calculation of the local exponent H. So by analogy with conventional definition of the local temperature in 
statistical physics one can generalize the conception of the Hurst exponent and use it for short samplings. The reasons 
for changing of the Hurst exponents can be the following: a) slow changing of type of correlations inside the noises; 
b) the presence of the regular signal inside the noises. So, in concrete applications the local Hurst exponent can serve 
as a quantitative characteristic reflecting the fractal properties of the EQ or TE event. It is obvious that the usage 
of long intervals (1000 registered points and more) for the calculation of the local Hurst exponents becomes useless 
and the important question is choosing the acceptable interval for calculating this parameter with high accuracy. The 
usual R/S analysis does not give the acceptable accuracy for the local Hurst exponent related to short samplings 
containing 100 — 120 points. So it is necessary to change the method of calculation of the local Hurst exponent for 
short samplings. The reliable calculation of the Hurst exponent averaged over short samplings turned out to be a 
nontrivial procedure and required elaboration of stable algorithms adjusted for averaging of short segments of the 
given samplings. 



We used another definition for the Hurst exponent |21j] , which turned out to be more effective for short samplings. 
The best results have been achieved in the usage of expression for the normalized dispersion, which relates differences 
of a random function to retardation time r: 

v( x ((B H (t + r)-B H (t)) 2 ) 2H 

V{T) = — wm — = r ' (27) 

here Bn{t) = J^., £j, is an integral random function. As a the result of numerical experiments it has became possible 
to calculate the local Hurst exponent with acceptable accuracy (4 — 5 decimal points) for samplings containing about 
80 — 120 registered points. 



IV. NON-MARKOV DISCRETE ANALYSIS OF SEISMIC DATA 



Here we will apply the discrete non-Markov procedure developed in the previous Section 2 for the analysis of the 
real seismic data. The basic problems, which we are trying to solve in this analysis, are the following. The first 
problem relates to a possibility of seismic activity description by statistical parameters and functions of non-Markov 
nature. The second problem relates to distinctive parameters and functions for differentiation of weak EQ's (with 
small magnitudes) from TE's. The third problem is the most important one and relates to strong EQ's forecasting. 
With this aim in mind we analyzed three parts of the real seismogram: before the event (EQ and TE), during the 
event and after the event. A typical seismogram contains 4000 registered points. The complete analysis includes the 
following information: phase portraits of junior dynamical variables, power spectra of four junior memory functions 
and three first points of statistical spectrum of non-Markovity parameter. We took into account also the values of 
numerical parameters characterizing the seismic activity. To analyze time functions we used also the power spectra 
obtained by the fast Fourier transform. The complete analysis exhibits great variety of data. 

We used 4 types of available experimental data courteously given by Laboratory of Geophysics and Seismology 
(Amman, Jordan) for the following seismic phenomena: strong EQ in Turkey (I) (summer 1999), a weak local EQ 
in Jordan (II) (summer 1998). As a TE we had the local underground explosion (III). The case (IV) corresponds to 
the calm state of the Earth before the explosion. All data correspond to transverse seismic displacements. The real 
temporal step of digitization r between registered points of seismic activity has the following values, viz., r = 0, 02s 
for the case I, and r = 0, 01s for the cases II-IV. The graphical information is classified as follows: 

Figs. 1-6 are referred to the case I; Figs. 7,8 correspond to the cases II and III considered together; Figs. 9, 10 illustrate 
the case IV. At first we consider the figures, which were obtained from the recordings corresponding to the states 
defined as before and during strong earthquake (EQ). Figs.l(a,b,c,d) (before EQ, state lb) and Figs.l(e,f,g,h) (during 
EQ, state I) demonstrate the temporal dynamics of four variables Wo(t), Wi(t), W2(t), Ws(t), which were calculated 
in accordance with Eqns. (|l4|), (|l6|). Let us note, that for convenience we use throughout initial variable Wo(t) as 
a dimensionless variable. From these figures it follows that for variable Wo the scale difference achieves the value 
more than 2500 (compare Figs. la and le)! In comparison with the cases (lb,c) the figures (lf,g) reveal the long-range 
and low frequency oscillations for variables W\ and Wi- One can calculate phase portraits in four-dimensional space 
of the obtained four dynamical variables Wo , W\ , W% , W3 as well. Figs. 2 and 3 show six projections on various 
two-dimensional planes of states: before (Figs. 2a,b,c, Figs. 3a,b,c) and during (Figs. 2d,e,f, Figs.3d,e,f) EQ. 

The phase portraits of the system analyzed demonstrate strong variations. The last arise owing to the transformation 
of the strained state of the Earth before the EQ to the state during the EQ. The most dramatic changes emerge in 
the phase plane (W 1: W ) (see Figs. 2a and 2d), plane (W 2 , W ) (Figs. 2b and 2e), plane (W 2 , W x ) (Figs. 3b and 3e) 
and(W3, W\) (Figs. 3a and 3d). One can notice strong qualitative variations in the structure of phase portraits in 
the following planes: (Wi, Wo), (W2, Wo) and (W2, W\). Besides, we can see the quantitative change of space scales 
of dynamic orthogonal variables. The plane projection (Wo,Wi) remind a strange attractor. The changes of phase 
portraits in other planes are less noticeable (compare Figs. 2c and 2f, Figs 3b, 3c, 3e and 3f). The weakest change 
is revealed in the phase portrait in the plane (W3, W2). Probably, this phase portrait is less informative and encloses 
quasi-invariant part of the total phase portrait. Besides the spatial scales change of the orthogonal variables W3 and 
W2 , other essential deformations of this phase portrait were not observed. 

As it has been mentioned above it is convenient to analyze the power spectra for comparison of memory functions. 
One can divide these spectra into the following regions: ultra-low frequency range (ULFR), low frequency range (LFR), 
middle frequency range (MFR) and high frequency range (HFR). Figs. 4 and 5 demonstrate spectra of four memory 
functions Mo, Mi, M2, M3 before and during EQ. Before (Fig. 4a) and during EQ (Fig. 4c) the power spectrum of 
the initial TCF Mo has a fractal form l/uj a in double- log scale. One can observe a peak in ULFR (Fig. 4c) during 
EQ. The power spectra of the first and second memory functions during EQ (Figs. 4d and 5c) have also the fractal 
structure. The last one reflects the existence of linear frequency dependence in double- log scale within the LFR, MFR 



and HFR. The similar fractal-like behavior for the Turkish strong EQ is preserved for the third memory function for 
the state during EQ (see Fig. 5d). 

Figs. 6 demonstrates the power spectra of the first three points of the statistical spectrum of non-Markovity 
parameter for the states before (a, b, c) and during (e, f, g) strong the EQ. One can make the following conclusions 
from Figs. 6 a-d. On the first level of relaxation process (see, Fig. 6a) the strained state of the Earth crust before 
EQ can be associated with Markov and quasi-Markov behavior in ULFR and LFR, correspondingly. The influence 
of non-Markov effects is reinforced in MFR with 5 • 10^ 2 f.u. < u> < f.u., (1 f.u. = 2tt/t). Strong non- Markovity 
of the processes considered for £i(w) takes place in HFR with lQ^ 1 f.u. < w < 0.5/.U.. Simultaneously we have the 
numerical values £2(0;), £3(0;) ~ 1 in the whole frequency region(see, Figs. 6 b, c). But this behavior implies that 
strong non-Markovity effects are observed in these cases. 

The similar picture becomes unrecognizable for seismic state during the strong EQ (see, Figs. 6 d-f). Firstly, it is 
immediately obvious that £i(w) ~ 1 on first relaxation level. 

Secondly, the second and third relaxation levels are non-Markovian (see, Figs. 6 c, f). Thus, the behavior of seismic 
signals during the strong EQ is characterized by strong non-Markovity on the whole frequency region. 

Figs. 7 depicts power spectra of MF Mq and Mi for seismic states II and III. Figs. 8 shows power spectra of the first 
three points of non-Markovity parameter Si(u>),i = 1,2,3. The preliminary results suggest that there is remarkable 
difference between weak EQ's and TE's especially in the area of low frequencies. 

The analysis of the phase portrait for weak EQ's and underground TE's lead to the following conclusions. Firstly, 
these portraits cannot be differentiated. It can be seen from the range of spatial scales of the dynamical variables 
Wi and Wj and from the analysis of the phase points distribution forms. Secondly, it is necessary to remark some 
peculiarities in power spectra of /Uj(w), i = 0,1,2,3 (see, Figs. 8 ) for the cases II and III. All these spectra have 
distinctive similarities for the memory functions M,(i) with numbers i = 0, 1, 2 and 3. The character and the form of 
the spectra considered for the cases II and III are very similar to each other. The same similarity is observed for the 
three non-Markovity parameters ei(u)), i = 1,2,3 (see, Figs. 8). 

Nevertheless the analysis of the power frequency spectra allows to extract distinctive specific features between the 
weak EQ's and the TE's. Such quantitative criteria can be associated with frequency spectra of memory functions 
Hi(u>) characterizing the long-range memory effects in seismic activity. This new criterion allows to tell definitely a 
weak EQ from a TE, viz., to differentiate case II from case III. 

A close examination of Figs. 8a and 8d shows that this distinction appear in frequency behavior of the first point 
of non-Markovity parameter £i(w) close to the zero point oj = 0. Specifically, the ratio of values £i(0) for weak EQ 
and TE equal £{ 7 (0)/£{ /7 (0) = 0.92/0.57= 1.61. 

Let us to analyze the results of seismic activity characterizing the calm Earth state. Figs. 9,10 present the results 
of this analysis. They will be useful for the comparison with the results obtained for EQ's and TE's. The projections 
of the phase cloud on all six planes (Wi,Wj), i ^ j exhibit approximately the similar distribution of phase points. 
The power spectra for the memory functions with the same parity (see, Figs. 9) have a similar form. For example, 
for even order functions /Uo(w) and ^(w) one can notice sharp peaks near the frequency 0.2 f.u. (see, Figs. 9a, c). 
In the spectrum of the senior function /^(w) (see Fig. 9c) additional peaks in HFR appears. One can notice two 
groups of characteristic peaks near 0.2 f.u. and OAf.u in odd memory functions /ii(w) and ^3(0;) (Figs. 9b, d). With 
the increase of order of the memory function the pumping over effect of peak intensities from the MFR to the HFR 
takes place. The frequency behavior of the three points of non-Markov parameters ei(w), 62 (u>) and €3(0;) appeared 
to be practically the same. The behavior of the functions Ci(u>) exhibits the typical non-Markov character with small 
oscillations of random nature at LFR. The spectral characteristics of the system IV are very useful in comparison to 
the results obtained for the system I (before strong EQ). 

Our observation shows that the zero point values of non-Markovity parameters for calm Earth state are equal 
s{ v (0) : £^(0) : £3^(0) « 4.99 : 0.947 : 0.861. These values are convenient for the comparison with similar values 
for the Earth seismic state before the strong EQ: £{(0) : 4(0) : 4(0) w 214.3 : 0.624 : 0.727. The change of ratio of 
the two first non-Markovity parameters £i(0)/£2(0) is particularly striking . This ratio is equal to 5.27 for the calm 
Earth state, then it comes into particular prominence for the state before strong EQ: £{ (0)74(0) w 343.4. Thus, this 
ratio changes approximately in 60 times! Hence, the behavior of this numerical parameter is operable as a reliable 
diagnostic tool for the strong EQ prediction. The foregoing proves that the indicated value drastically increases in 
process of nearing to strong EQ. 

Finishing this Section, we give some preliminary suggestions relating to the strong EQs forecasting. They are related 
in comparison of frequency spectra obtained for the calm Earth (Figs. 11,12) and seismic activity data registered 
before a strong EQ (see, Figs. 2a,b,c; 3a,b,c; 4a, b; 5a, b and 6a, b). The comparison of the phase portraits demonstrates 
the following peculiarities. In the phase portraits calculated for the senior dynamical variables (W2, Wi), (W3, Wi) 
and (Ws,W2) obtained for cases I and IV the distinctions arc not noticeable (see Figs. 3a, b). These distinctions 
become noticeable in the phase portraits of junior variables (Wi,Wo), (W2,Wq) and (W3,Wo) (see, Figs. 2a,b,c). 
One can observe a gradual stratification of the phase clouds with the growth of elastic deformations before the strong 



EQ. It is necessary to recall the double frequency difference for systems I and IV when comparing the frequency 
plots. The dependence ^o( w ), Mi(^) 5 M2(w) and ^3(0;) for systems I and IV (see Figs. 4a, b; 5a, b; and 14 a-d) is 
approximately similar, and qualitative difference is not noticeable. One can notice some visual difference only for two 
spectra: for the third memory function spectrum Hs(u)) and for the ULFR of the memory function fj,o(u>). So the 
power spectra of memory functions can be used for the strong EQ forecasting. One can notice the similar changes 
in the behavior of the functions ei(o>), e 2 (o>) and e 3 (o>) (see, Figs. 6 a,b,c and 10 a,b,c). So one can conclude that 
careful investigations of frequency behavior of memory functions Hi(u>) and functions ei(u>) describing the statistical 
non-Markovity parameters provide an accurate quantitative method of the strong EQ's forecasting. It is necessary 
to investigate carefully the power spectra with the accurate localization of an object and source, generating seismic 
signals, for further elaboration of this method. 

For more complete understanding of non-Markov properties of seismic signals we give some kinetic parameters of 
our theory in Tables I-III. In Table I the full sets of kinetic parameters describing non-Markov stochastic processes 
in five various seismic states have been presented: before strong EQ(Ib), during strong EQ (I), during weak EQ (II), 
during TE (III) and for the calm Earth state (IV). The data cited in this Table are indicative of non-equilibrium 
properties (parameters Ai, A2 and A3), long-range memory effects (parameters Ai and A2) and non-Markov peculiarities 
(parameters £i(0), £2(0) and £3(0)). The differences under observation for various seismic states are sufficient to allow 
definite conclusions. 

For purposes of clarity, Table II illustrates the of comparison of specific kinetic non-Markov parameters for two 
seismic state: before strong EQ(Ib) and calm Earth states (IV). As will be seen from Table II differences of parameters 
for these two states vary within a broad range: from 2.8 (parameter A 2 ) to 44.0 (for parameter e 1 (0)). Similarly, 
Table III contains comparison data for the other two seismic states: during weak EQ(II) and during underground 
TE(IV). Differences of parameters in this case are established within more narrow limits: from 2.486 (for parameter 
Ai) to 1.614 (for parameter £i(0)). 

Thus, the existence of discreteness and long-range memory in behavior of seismic signals opens up new fields of use 
in the analysis of the Earth seismic activity. We can state with assurance that the differences under observation favor 
the view that the non-Markov parameters of our theory will be available for strong EQ's forecasting and differentiation 
of TE's from weak EQ's. 

V. LOCAL HURST EXPONENT CALCULATIONS FOR AVAILABLE EARTHQUAKES AND 

TECHNOGENIC EXPLOSIONS DATA 

Available data for the calculation of the local Hurst exponents contains 3000-5000 registered points describing the 
visible part of a wavelet. This number of the recorded points allows to use the procedure of the local Hurst exponent 
H(t) calculation. For the realization of the procedure described in the previous Section 3 it is necessary to divide 
the whole sampling containing 25000 points into small intervals of 100-200 points, where the local Hurst exponent is 
supposed to be constant. On Fig. 11 we show a typical plot of the function H(t) calculated for a typical EQ. The 
same features of H(t) behavior are conserved for wide class of available weak EQ's. Then we obtained the calculated 
values of the local Hurst exponents H{t) for available EQ's and explosions. Figs. 12 exhibits the typical behavior of 
these functions. The sharp decreasing (0.1) of the local Hurst function during "an event" is typical for explosions. 
Then the values of the function H{t) are relaxing slowly to their initial values. For EQ's one can notice more gradual 
change of H(t) before the event. The relaxation of H{t) starts from higher (0.85) values and it comes back faster 
to its initial values in comparison with explosions. Such behavior is preserved for weak signals, when the ratio S/N 
decreases. For these cases the criterion of EQ or TE distinction is related to the amplitude of the Hurst exponent 
change during the analyzed event. It is necessary to increase the number of registered points per unit of real time in 
order to obtain a more distinctive picture, which can be more useful in differentiation of these events. It is related to 
the fact that the sensitivity of correlations of a random fractal value changing is associated with the lower temporal 
limit of the corresponding measurements. The smoothed change of H{t) obtained for EQ's opens a possibility of more 
accurate registration of H{t) before the visual wavelet of EQ's. 

VI. CONCLUSION 

We want to stress here again that these presented method have been applied successfully for differentiation of EQ's 
from TE's. We hope that the results of this analysis can be applied to a set of phenomena related with differentiation 
of similar signals of different nature. In the result of this analysis we received a new possibility to forecast of strong 



EQ's approaching, analyzing only seismograms recorded for transverse seismic waves. Secondly, we received sufficient 
amount of information for the definite differentiation of weak EQ's from TE's. 

In this paper we have presented the results of application of two new methods for the study of dynamic, kinetic and 
spectral properties of seismic signals depicting EQ's and TE's modulation. By the discrete non-Markov stochastic 
processes and the local Hurst exponent analysis we have found explicitly some features of several different states of 
the Earth crust: states of the the Earth before and during strong and weak EQ's, during TE's. The used methods 
allow to present the seismogram analyzed in the form of set non-Markov variables and parameters. They contain 
great amount of the qualitative and quantitative information about seismic activity. 

The dynamic information is contained in time recordings of new orthogonal dynamic variables, different plane 
projections of the multidimensional phase portrait and the time dependence of local Hurst exponent. The information 
on the kinetic, spectral and statistical properties of the system is expressed through time dependence of the initial 
TCF, memory functions of junior orders, their power and frequency spectra of the first three points of statistical 
spectrum of non-Markovity parameter. 

The main advantage of our two new methods is a great amount of supplementary information about the properties 
of seismic signals makes. The problem is its correct application. What kind of possibilities can one expect? It is 
possible to answer as follows. Firstly, our preliminary study, convincingly demonstrates that the relevant and valuable 
information on non-Markov and discrete properties of the system considered is contained in seismic signals. In all the 
studied systems (I- IV) we have found out unique manifestations of Markov, quasi-Markov and non-Markov processes 
on the particular behavior of the signals in a broad range of frequencies. 

The similar results cannot be obtained, on principle, by other methods used in the analysis of seismic activity. 

Secondly, in the non-linear non-Markov characteristics some of well-known spectral effects are evident. Among them 
the following effects arc exhibited noticeably: fractal spectra with an exponential function cj~ a 7 which are connected 
to the phenomenon of usual (SOC) and restricted (RSOC) self-organized criticality [15,16], behavior of some frequency 
spectra in the form of white and color noises. Thirdly, the frequency spectra introduced above are characterized by 
the particular alternation of Markov (fractal) and non-Markov spectra (such as color or white noises). The similar 
alternation resembles in particular the peculiar alternation of effects of a non-Markov and non-Markov behavior for 
hydrodynamic systems in the statistical physics of condensed matter detected in papers [17], [18] for the first time. 
The fine specification of such alternation appears essentially different for studied states I- IV. These features allow to 
view optimistically the solution of the problem of forecasting strong EQ's and differentiation TE's from weak EQ's. 
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VIII. FIGURE CAPTIONS 



Figure 1. The temporal dynamics of the first four dynamic variables Wo(t), W\(t), W2(t), W 3 (t): a,b,c,d - before 
strong EQ, e,f,g,h - during of strong EQ. During strong EQ fluctuation scale increase drastically. It makes up: 2.5 *10 3 
for initial variable Wo(i),10 2 for the first orthogonal variable Wi(t),10 for W2(t) and 2 for W 3 (t). The existent trend 
vanishes gradually at transition from the initial variable W (t) to the third orthogonal variable W 3 (t) . The fluctuation 
scale decreases sharply during the strong EQ. 

Figure 2. The phase portrait projections on the planes of orthogonal variables W ,Wi (a), W ,W 2 (b), W , W 3 
(c)- before the strong EQ(Ib) and W , Wi (d),W , W 2 (e), W a ,W 3 (f) during the strong EQ (I). The sharp difference 
is distinct for seismic states lb and I. The randomization of the phase portrait for state I begins from plane Wq, W 2 . 
Together with the difference of the scale of fluctuation, one can observe the asymmetric distribution of phase clouds 
everywhere. 

Figure 3. The phase portrait projections on the planes of orthogonal variables Wi,W 2 (a), Wi, W 3 (b), W 2 ,W 3 
(c)- before the strong EQ (lb) and on planes Wi,W 2 (d),Wi, W 3 (e), W 2 ,W 3 (f) during of strong EQ (I). All phase 
clouds for seismic state lb are symmetrical as opposed to Figs. 2. Sharply marked asymmetry and stratification of 
phase clouds, what resembles known situation for myocardial infarction in cardiology, are observed for state I (d,e, 
and f). 

Figure 4. The power spectra of the two first memory functions /io and fi\ : (a,b) -before the strong EQ(Ib), (en- 
during of the strong EQ (I). For the cases a,c and d we observe fractality and self-organized criticality (SOC). SOC 
exists for the whole frequency range for state lb. However, we observe restricted SOC at c and d cases only in 
frequency range down to 2.5 * 10~ 3 units of (2-k/t). Restricted SOC is characterized by sharp decreasing of intensity 
on two orders for c and d cases! One can see color noises nearby 0.1 and 0.2 f.u. for in state lb. 

Figure 5. The spectra of two memory functions \x 2 and \x 3 : (a,b) -before the strong EQ, (c,d)- and during the strong 
EQ. One can observe color noises in cases a,b and d. Fractal-like spectrum on ultra-low frequencies is appreciable in 
addition to cases c and d. The spectra for states lb and I are sharply different from each other both to intensity and 
to spectral peaks positioning. 

Figure 6. Frequency spectra of the first three points of non-Markovity parameters e\ e 2l e 3 : (a,b,c)- before the strong 
EQ, (d,e,f)- during the strong EQ. Markov and quasi-Markov behavior of seismic signals is observed only for e\ in 
state lb. All remaining cases (b,c,d and d) relate to non-Markov processes. Strong non-Markovity is typical for cases 
b , c (state lb) and for case d (state I). In behavior of e 2 (oj) and e 3 (td) one can see a transition from quasi-Markovity 
(at low frequencies) to strong non-Markovity (at high frequencies). 

Figure 7. The power spectra for the two first memory functions fj,o and fii: (a,b) -during weak EQ, (c,d)- during 
technogenic explosion. In cases lb and I the spectra are characterized by strong differences especially on ultra low 
frequencies. They have very low intensity for /zo on low frequencies (cases a and c) and color like behavior for Hi for 
states II and III(cases b and d). Unexpected peaks exist in system III in LFR. Color and intensity distribution of the 
spectra is different for states II and III. 

Figure 8. The frequency spectra of the first three points in statistical spectrum of non-Markovity parameters ei, e 2l 
e 3 : (a,b,c)- during weak EQ, (d,e,f)- during TE . All spectra are characterized by strong expressed non- Markovity ( 
€i ~ 1) for the whole frequency range. Weak quasi- Markovity is observed near zero frequency for cases a and d (ei 
vary from 0.5 up to 6.5). A noticeable difference for states II and III exists in behavior ei(uj) in point lo = 0. Due to 
this fact, one can develop reliable approach to differentiation between weak EQ's and underground TE's. 

Figure 9. The power spectra of memory functions /io(^), M2(w) an d M3( w ) f° r the calm state of the Earth 

before explosion. All functions Hi(u>), i = 0,1,2,3 have approximately similar fractal behavior with restricted SOC 
and color noises close to 0.2 and 0.4 f.u. The maximum of intensity emerges close to the frequency 4 * 10~3/.u. . A 
slight change and redistribution of intensity of power spectra occur with the increase of order of memory function. 

Figure 10. The power spectra of the first three points in statistical spectrum of non-Markovity parameter ei, e 2 , e 3 
for calm state of the Earth before explosion (IV). Due to similar frequency behavior of all memory functions the 
functions £j(w) , i=l,2 and 3 have approximately similar frequency behavior and therefore demonstrate strong non- 
Markovity on all levels. The initial parameter €\{lo) is non-Markovian with the exception of slight quasi-Markovity 
close to low frequencies below 0.1 f.u. As a result of this the possibility appears for forecasting the strong EQ's by 
registration of disappearance of strong non-Markovity and appearance of pronounced Markov time effects. 

Figure 11. The typical temporal behavior of the Hurst exponent H(t) calculated for EQ's. One can see sharp 
decreasing of H(t) on 15 % during EQ. After that a gradual restore of the Hurst exponent H{t) to normal value « 1 
takes place. 

Figure 12. The comparative analysis of the Hurst exponent H{t) behavior during the weak EQ (a,c) and for the TE 
(b,d). During the weak EQ's one can see sharp decreasing of H(t) on 15 % and almost 90 % during the TE. These 
observations enable us to develop a new approach to differentiation the TE's from weak EQ's. 



Table I. Set of kinetic non-Markov parameters of discrete stochastic processes in various seismic states 





before strong 
EQ, lb 


during strong 
EQ, I 


during weak 
EQ, II 


during TE, III 


calm state of 
Earth, IV 


Ai, units of r _i 


-0.0052275 


-0.00010709 


-0.32465 


-0.17203 


-0.22972 


A2, units of r -1 
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Ai, units of 
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0.059232 


0.021134 


A2, units of t~ 2 


0.31541 


4.5948e-005 
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-0.032079 


0.11266 


6,(0) 
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Tabic II. Comparison of kinetic non-Markov parameters for two seismic states: before strong EQ (lb) and calm 

Earth state (IV) 



ratio of 
parameters 


e{ b (0)/s{ v (0) 




A{ b /A{ v 




numerical value 


42.94 


1:22.0 


1:1.3 


1:0.7 
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Table III. Comparison of kinetic non-Markov parameters for two seismic states: during weak EQ (II) and during 

underground TE (III) 



ratio of 
parameters 






eV(0)/e{"(0) 


numerical value 


2.486 


1.887 


1.614 
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